$ontext
This model uses the results of regional models of with the same scenario name to
find the supply curve of biojet fuel from miscanthus or switchgrass planted on the
top 25% of marginal lands. The model assumes a $50/Mg CO2-e carbon price.
$offtext

$set scenario E8_M025_50C_SC

Option   Solprint = off;
Option   Limrow = 0;
Option   Limcol = 0;


sets      plot   biomass source locations  /
$include wsc_pseudocounties_E8_025.csv
/
         brfn    biorefinery sites /
$include brfn_all.csv
/
         feed    /miscanthus, switchgrass/
         crop /ncrp, soy, srg, wht, pnt, ric, osg, cot, hay, oat, crn, brl, bns, sbt, sfl/
         crop_cat /ncrp, crp/
         extra   feedstock classification        /extra1, extra2/
         comp    feedstock composition component /cellulose, hemicellulose, lignin, HHV/
         tech         tech type   /lce, ft_jet/
         term            terminal locations  /
$include terminal_list.csv
/
         trans        transportation mode /hours, dist, total_cost/;

set      state           /
$include state_list.csv
/;


alias (tech, tech2);

* Regional links file is aggregated from regional model results with the cost of
* transport between pseudocounties and potential biorefineries
Parameters       pbcost(plot, brfn)        '$ per ton for feedstock transport'  /
$ondelim
$include regional_links_E8_MS025.csv
$offdelim
/

                 btcost(brfn, term)  'cost of transport of fuel to end-use'

         acres(plot, crop) 'Land area available' /
$ondelim
$include wsc_acres_E8_025.csv
$offdelim
/
         wta(plot, feed, crop)    'Willingness to accept prices for biomass'    /
$ondelim
$include wsc_wta_E8_025.csv
$include wsc_wta_E8_025_sg.csv
$offdelim
/
         yield(plot, feed)  'yields of energy crops'     /
$ondelim
$include wsc_yields_E8_025.csv
$include wsc_yields_E8_025_sg.csv
$offdelim
/

         feed_MC(feed)                   moisture content of dry biomass feedstocks   /
miscanthus           0.15
switchgrass              0.15/

         harvest_loss            fraction of total biomass yield not harvested    /0.1/
         logistics_loss          fraction of biomass lost in logistics           /0.05/

         conv_eff(feed, tech)            conversion efficiency for converting feedstock to fuel

         conv_factor(tech, comp)    converts conversion efficiencies into gallons product per ton feedstock             /
lce.hemicellulose    176.9
lce.cellulose        172.85
ft_jet.HHV        7.102/
         tech_eff(tech, comp)       conversion efficiency for cellulosic technologies     /
lce.hemicellulose    0.765
lce.cellulose        0.72
ft_jet.HHV        0.465/

         electricity(tech, feed)               electricity produced per ton of biomass
         elec_eff(tech)          /
lce              0.0366
ft_jet        0.0032/
         naptha(tech, feed)                    naptha produced per ton of biomass
         naptha_eff(tech)                /
ft_jet        0.0/
         gge_conversion(tech)          converts fuel quanitities into equivalent gallons of gasoline   /
lce              0.657
ft_jet        1/

         vmt_fraction(term)      'fraction of 2015 vmt allocated to each terminal'  /
$ondelim
$include terminal_vmt_fraction.csv
$offdelim
/
         fuel_price_base              fuel price for each model run                                  /2.8/
         fuel_price(term)        terminal specific fuel price accounting for historic price differentials
         naptha_price            price paid for co-product naptha
         elec_price              price of electricity in dollars per kWh                         /0.05/
         carbon_price            price of carbon in dollars per metric ton CO2e                  /50/
         total_supply            target supply quantity for cost min model                      /3000/
         marg_cost_land(crop)    adder to the cost of land conversion as a linear function of converted land
         distance(plot,brfn)
;

naptha_price = 0.86*fuel_price_base;
*adder for additional land rent on cropland which is $5 per thousand acres
marg_cost_land(crop) = 5/1000/10000;
marg_cost_land('ncrp') = 0;
yield(plot, feed) = yield(plot,feed)*(1 - harvest_loss);
wta(plot,feed,crop) = wta(plot,feed,crop)/(1 - harvest_loss);

Table    biomass_composition(feed, comp)
$ondelim
$include feedstock_composition.csv
$offdelim
;

Table    feedstock_od_table(plot, brfn, trans)
$ondelim
$include p2b_limited_025.csv
$offdelim
;

Table        terminal_od_table(brfn, term, trans)
$ondelim
$include wsc_brfn2term.csv
$offdelim
;
Table            terminal_cost(term, tech)
$ondelim
$include terminal_cost.csv
$offdelim
;

Parameter state_terminals(term, state)      /
$ondelim
$include terminal_fips.csv
$offdelim
/
         state_brfn(brfn, state)      /
$ondelim
$include brfn_fips.csv
$offdelim
/
         state_price_delta(state)        /
$ondelim
$include state_price_delta.csv
$offdelim
/
         crp_cat_matching(crop_cat, crop)        /
$ondelim
$include crop_cat_crosswalk.csv
$offdelim
/
         regional_locations(brfn, tech, extra)   /
$ondelim
$include regional_locations_MS.csv
$offdelim
/
;

btcost(brfn, term) =  terminal_od_table(brfn, term, 'total_cost');
distance(plot,brfn) = feedstock_OD_table(plot,brfn,'dist');

fuel_price(term) = fuel_price_base;
* To make the fuel prices spatially disaggregated, you can and a state specific price delta with the following line.
* + sum(state, state_terminals(term, state)*state_price_delta(state));



conv_eff(feed, tech) = sum(comp, biomass_composition(feed, comp)*tech_eff(tech, comp)*conv_factor(tech, comp))/1000;
electricity(tech, feed) = elec_eff(tech)*biomass_composition(feed, 'HHV')/1.1023/.0036/1000;
naptha(tech, feed) = naptha_eff(tech)*biomass_composition(feed, 'HHV')/1.1023/0.1216089/1000;

*Carbon accounting parameters
Parameters
* LUC carbon from GREET
         LUC_CO2(plot, feed, crop)
         LUC_by_cat(plot, feed, crop_cat)    'LUC emissions from GREET Mg C/acre'    /
$ondelim
$include C_dLUC.csv
$offdelim
/
         cultivation_co2(plot, feed)       'biomass cultivation emissions kg CO2e/dry ton sold'
         c02_per_ton_mile        'Carbon emissions in biomass transport Mg CO2e/dry ton-mile'  /0.000177/
         conversion_co2(tech)    'Carbon emissions at brfn Mg CO2e/dry ton consumed'    /
lce              0
ft_jet        0.01775/
         CO2_elect(brfn, tech)   'Carbon emissions of displaced electricity kg/MWh'
         state_elect_CI(state)   'Carbon intensity of electricity in each state in kg/MWh from 2019 eGRID' /
$ondelim
$include state_electricity_CI.csv
$offdelim
/
         fuel_transport_CO2      'Carbon emissions kg per 1,000 gallon-miles of fuel transport' /0.645/
         bt_CO2(brfn,term)       'fuel transport emissions kg CO2e/1,000 gal-miles'
         end_use_CO2(tech)       'Carbon emissions at the end use kg/1,000 gallons'  /
lce              0.0366
ft_jet        0.0032/
;

* Calculate emissions over each fuel link
bt_CO2(brfn, term) =  fuel_transport_CO2*terminal_od_table(brfn, term, 'dist');

* Assign state electricity CI to brfns in the state, convert to Mg/Mwh
CO2_elect(brfn, tech) = (sum(state, state_brfn(brfn, state)*state_elect_CI(state)))/1000;

* Calculate cultivation emissions based on yield from crop budgets and GREET1_2021, convert to Mg/ton
cultivation_co2(plot, 'switchgrass')$(yield(plot, 'switchgrass') gt 0.01) = (38.56 + 13.56/yield(plot, 'switchgrass'))/1000;
cultivation_co2(plot, 'miscanthus')$(yield(plot, 'miscanthus') gt 0.01) = (28.37 + 2.97/yield(plot, 'miscanthus'))/1000;

*Convert soil organic carbon changes into CO2 emissions
LUC_CO2(plot, feed, crop) = - sum(crop_cat, crp_cat_matching(crop_cat, crop)*LUC_by_cat(plot, feed, crop_cat))*44/12;


*Removing unneeded varibales and equations
sets spar(plot, feed, brfn)  set used to eliminate unnecessary links
     spar1(plot, feed, crop)       set used to eliminate unnecessary supplies
         spar2(brfn, tech, extra)
         spar4(feed, brfn, tech, extra)
         spar3(brfn, tech)
         spar5(brfn, term, tech)
         spar6(feed, brfn)
         spar7(plot, crop)
         spar8(plot, feed);

spar(plot, feed, brfn)$(pbcost(plot, brfn)*sum(crop, acres(plot, crop))*yield(plot, feed) > 0.5)=yes;
spar1(plot, feed,  crop)$(acres(plot, crop)*sum(brfn, pbcost(plot, brfn))*yield(plot, feed) > 0.5)=yes;
spar2(brfn, tech, extra)$(regional_locations(brfn, tech, extra) gt 0.5) = yes;
spar3(brfn, tech)$(sum(extra, regional_locations(brfn, tech, extra)) gt 0.5) = yes;
spar4(feed, brfn, tech, extra)$(regional_locations(brfn, tech, extra) gt 0.5) = yes;
spar5(brfn, term, tech)$(sum(extra, regional_locations(brfn, tech, extra)) gt 0.5) = yes;
spar6(feed, brfn)$(sum((tech, extra), regional_locations(brfn, tech, extra)) gt 0.5) = yes;
spar7(plot, crop)$(acres(plot, crop)*sum(brfn, pbcost(plot, brfn))*sum(feed, yield(plot, feed)) gt 0.5)=yes;
spar8(plot, feed)$(sum(crop, acres(plot, crop))*sum(brfn, pbcost(plot, brfn))*yield(plot, feed) > 0.5)=yes;


*Assign missing links a very high cost
pbcost(plot, brfn)$(pbcost(plot, brfn) lt 0.05) = 200;


*Fix the parameters to have units of $10M, 10M gallons, 10k tons and 10k Mg Co2e
pbcost(plot, brfn) = pbcost(plot, brfn)/1000;
acres(plot, crop) = acres(plot, crop)/10000;
btcost(brfn, term) = btcost(brfn, term)/100;
wta(plot, feed, crop) = wta(plot, feed, crop)/1000;
carbon_price = carbon_price/1000;

Positive Variables
         adopted_acres(plot, feed, crop)         area planted in energy crop "feed" in "plot" with incumbent "crop"
         biomass_consumed(feed, plot, crop)      annual quantity of feedstock consumed from plot "i" at price level "plev"
         p2b(feed, plot, brfn)                   annual quantity of biomass from plot "i" taken to biorefinery "k"
         Qb(brfn, tech, extra)                   quantity of fuel "b" produced at biofinery "k"
         T(brfn, term, tech)                     quantity of fuel delivered to terminal from biorefinery
         x1(feed, brfn, tech, extra)             linearizing variables for biorefinery costs
         crop_displaced(crop)                    total displaced area for a given crop
*         total_supply
;

Variable
         total_cost
         total_carbon
         profit;

binary variable
         xi(brfn, tech, extra);

Equations
         annual_cost                     total annual cost for biofuel supply
         annual_profit
         xi_constraint(brfn, tech, extra)       constraint for integer variable xi
         Qb_constraint(brfn, tech, extra)       relates Qb to x1's

*Constraints to set up biofuel networks
         land_constraint(plot, crop)             limits conversion to available land and prevents double counting of land
         supply_constraint(plot, feed, crop)     limits biomass leaving a plot at a price level to what is available
         plot_flow(feed, plot)                   constrains feedstock leaving plot to what is consumed
         brfn_cap(feed, brfn)                   limits fuel production capcity at biorefinery to be less than or equal to the feedstock inputs
         brfn_flow(brfn, tech)                   limits the fuel leaving the brfn to be less than the fuel produced
*         total_demand                            aggregates the fuel quantities to the set limit
         crop_accounting(crop)
         carbon_accounting

;


Parameters
         ac(tech)                         fixed cost component in linearized biorefinery costs      /
lce              3.6917481
ft_jet        33.1629/
         bc(tech)                         feed capacity depedent component in linearized biorefinery costs      /
lce              0.4528
ft_jet        0.48858/
         iir             /0.1/
         lifetime        /20/
         CRF
         ao(tech)                         fixed cost component in linearized biorefinery costs      /
lce              0.4356352
ft_jet        3.15/
         bo(tech)                         tech cap dependent component in linearized biorefinery costs    /
lce              0.05952
ft_jet        0.0764/
         MM(tech)                               Maximum capacity                                  /
lce              131
ft_jet        200/
         a(brfn, tech, extra)
         b(brfn, tech, extra)
         M(brfn, tech, extra);

CRF = iir*(1+iir)**lifetime/((1+iir)**lifetime - 1);


a(brfn, tech, extra) = (CRF*ac(tech) + ao(tech));
a(brfn, tech, 'extra2') = (CRF*ac(tech) + ao(tech)) + 0.1;
b(brfn, tech, extra) = (CRF*bc(tech) + bo(tech));
M(brfn, tech, extra) = MM(tech);

*Objective
annual_cost..            total_cost =e= sum((spar(plot, feed, brfn)), pbcost(plot, brfn)*p2b(feed, plot, brfn)) + sum((spar1(plot, feed, crop)), wta(plot, feed, crop)*biomass_consumed(feed, plot, crop))
                         + sum((feed, brfn, tech, extra), b(brfn, tech, extra)*x1(feed, brfn, tech, extra)) + sum((brfn, tech, extra), a(brfn, tech, extra)*xi(brfn, tech, extra))
                         + sum((brfn, term, tech), T(brfn, term, tech)*(btcost(brfn,term)+ terminal_cost(term, tech))) + sum(crop, marg_cost_land(crop)*crop_displaced(crop));

annual_profit..          profit*1000 =e= sum((brfn, term, tech), fuel_price(term)*gge_conversion(tech)*T(brfn, term, tech))+ elec_price*sum((brfn, tech, extra, feed), electricity(tech, feed)*x1(feed, brfn, tech, extra))
                         + naptha_price*sum((brfn, tech, extra, feed), naptha(tech, feed)*x1(feed, brfn, tech, extra)) - total_cost - carbon_price*total_carbon/1000;

*Subject to:
land_constraint(plot, crop)$spar7(plot, crop)..    sum(feed, adopted_acres(plot, feed, crop)) =l= acres(plot, crop);
supply_constraint(plot, feed, crop)$spar1(plot, feed, crop)..     biomass_consumed(feed, plot, crop) =e= yield(plot, feed)*adopted_acres(plot, feed, crop);
plot_flow(feed, plot)..        sum(spar1(plot, feed, crop), biomass_consumed(feed, plot, crop)) =g= sum(spar(plot, feed, brfn), p2b(feed, plot, brfn));
brfn_cap(feed, brfn)..        sum(spar(plot, feed, brfn), p2b(feed, plot, brfn)*(1 - logistics_loss))
                                                                 =g= sum((tech, extra), x1(feed, brfn, tech, extra));
brfn_flow(brfn, tech)..          sum(term, T(brfn, term, tech)) =l= sum(extra, Qb(brfn, tech, extra));
crop_accounting(crop)..            sum(spar8(plot, feed), adopted_acres(plot, feed, crop)) =e= crop_displaced(crop);

*Demand constraint for use in cost minimization models
*total_demand..           sum((brfn, term, tech), gge_conversion(tech)*T(brfn, term, tech)) =g= total_supply;

xi_constraint(brfn, tech, extra)..      sum(feed, x1(feed, brfn, tech, extra)) =l= M(brfn, tech, extra)*xi(brfn, tech, extra);
Qb_constraint(brfn, tech, extra)..      Qb(brfn, tech, extra) =l= sum(feed, x1(feed, brfn, tech, extra)*conv_eff(feed, tech));

carbon_accounting..      total_carbon =e= sum((spar1(plot, feed, crop)), LUC_CO2(plot, feed, crop)*adopted_acres(plot, feed, crop) + cultivation_co2(plot, feed)* biomass_consumed(feed, plot, crop))
                         + sum(spar(plot, feed, brfn), p2b(feed, plot, brfn)*distance(plot, brfn)*c02_per_ton_mile)
                         + sum(spar4(feed, brfn, tech, extra), conversion_co2(tech)*x1(feed, brfn, tech, extra)) - sum(spar4(feed, brfn, tech, extra), CO2_elect(brfn, tech)*electricity(tech, feed)*x1(feed, brfn, tech, extra))
                         + sum(spar5(brfn, term, tech), T(brfn, term, tech)*(bt_CO2(brfn,term)+ end_use_CO2(tech)));

*Limit biorefinery locations to those selected by the regional models
xi.up(brfn, tech, extra) = regional_locations(brfn, tech, extra);

adopted_acres.up(plot, feed, crop) = acres(plot, crop);

Option MIP = cplex;
Option iterlim = 1000000;
Option reslim = 18000;

Model USDA /ALL/;
*$ontext
FILE opt "Cplex option file" / cplex.opt /;
PUT opt;
PUT
'workmem 24000'/
'memoryemphasis 1'/
'brdir 1'/
'cuts    2'/
'probe   3'/
'parallelmode -1'/
'threads=0'/;
PUTCLOSE OPT;

USDA.optfile=1;
* Optimality criteria defines the allowable optimality gap as a fraction of the objective function.
USDA.OptCR = 0.01;
*$offtext
*Option savepoint=1;

*$ontext
Set      run  /run0*run20/;

alias  (extra,extras);
alias    (term, term2);

Parameter fprice(run)   defines the price for each run /
run0        3.8
run1        3.9
run2        4
run3        4.1
run4        4.2
run5        4.3
run6        4.4
run7        4.5
run8        4.6
run9        4.7
run10        4.8
run11        4.9
run12        5
run13        5.1
run14        5.2
run15        5.3
run16        5.4
run17        5.5
run18        5.75
run19        6
run20        7
/;


*$ontext
*Create results output files
file status_%scenario%;
file results_%scenario%   ;
file results_%scenario%_brfn;
*file results_%scenario%_feedstock_links;
*file results_%scenario%_fuel_links;
file results_%scenario%_pc;

*Make results files CSV
results_%scenario%.pc=5;
results_%scenario%_brfn.pc=5;
*results_%scenario%_feedstock_links.pc=5;
*results_%scenario%_pc.pc=5;
results_%scenario%_brfn.pw= 1500;
results_%scenario%.pw=1500;

$onend


$ontext
put results_%scenario%_feedstock_links;
put "scenario", "fuel_price ($/gge)", "source_id", "dest_id", "energy_crop", "quant_tons"/;

put results_%scenario%_fuel_links;
put "scenario", "fuel_price ($/gge)", "source_id", "dest_id", "pathway", "fuel_deliveries (MGY)"/;
$offtext

put results_%scenario%_pc;
put "scenario", "fuel_price ($/gal)", "source_id", "energy_crop", "avg_procurement", "quantity (bdt/yr)", "area (acres)", "LUC (Mg CO2e/yr)", "cultiv_GHG (Mg CO2e/yr)"/;
put results_%scenario%_brfn;
put "scenario", "fuel_price ($/gal)", "brfn_id", "technology", "extra", "fuel_output (MGY)", "electricity output (GWh/yr)", "naptha (MGY)", "feedstock_cap (bdt/day)",  ;
         loop (feed)     do
         put feed.tl,
         endloop;
put "capital_cost (M$)", "annual_cost (M$/yr)", "annual_capital (M$/yr)", "O&M (M$/yr)", "feedstock_procurement (M$/yr)", "feedstock_transport (M$/yr)", "fuel_distribution (M$/yr)",
         "carbon cost (M$/yr)", "ac ($/gal)", "max_trans_dist (miles)", "totalGHG (Mg/yr)", "CI (g/MJ)"/;

put results_%scenario%;
put "scenario", "fuel_price ($/gal)", "total_biofuel (MGal)", "total_land (Macres)", "total_GHG (mmMg/yr)";
         loop    (crop) do
         put     crop.tl,
         endloop;

         loop    (crop) do
         put     crop.tl,
         endloop;

         loop (tech) do
         put tech.tl,
         endloop;
put /;
$ontext
put " ", "$/gge", "MGGEY", "GWh/yr", "MGY",
         loop (tech, feed)$(feed_tech_match(feed, tech) gt 0.5)        do
         put feed.tl,
         endloop;
         loop (tech)     do
         put "# brfns",
         endloop;
         loop (feed)     do
         put "consumption ktons/yr",
         endloop;
put "billion $", "billion $",;
         loop (tech) do
         put "brfn capital (M$)",
         endloop;
put "million ton miles/yr", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr)", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr", "million ton miles/yr",
         "million acres", "milion acres"/;
$offtext

Parameter
                 feedstock_consumption(feed)
                 brfn_production(brfn, tech, extra)
                 brfn_consumption(brfn, tech, extra, feed)
*                 links(plot, brfn, feed)
*                 term_links(brfn, term, tech)
                 consumed_biomass(plot, feed, crop)
                 quant
                 biomass_production(crop)
*                 quant_last_run
                 brfn_quant(brfn, tech, extra)
                 brfn_count(tech)
                 elect_quant
                 naptha_quant
                 fuel_by_tech(tech)
                 total_annual_cost
                 brfn_capital_total(tech)
                 brfn_electricity(brfn, tech, extra)
                 brfn_naptha(brfn, tech, extra)
                 feedstock_cap(brfn, tech, extra)
                 capital_cost(brfn, tech, extra)
                 brfn_annual_cost(brfn, tech, extra)
                 annual_capital(brfn, tech, extra)
                 O_M(brfn, tech, extra)
                 feedstock_procurement(brfn, tech, extra)
                 avg_feed_pc(plot, feed)
                 feedstock_transport(brfn, tech, extra)
                 fuel_distribution(brfn, tech, extra)
*                 marginal_feedstock(brfn, tech, extra)
*                 mc(brfn, tech, extra)
                 avg_cost(brfn, tech, extra)
                 max_trans_dist(brfn, tech, extra)
                 crops_displaced(crop)
                 crop_acres_changed(plot, feed, crop)
                 total_land
                 total_GHG
                 brfn_GHG(brfn, tech, extra)
                 brfn_CI(brfn, tech, extra)
                 brfn_carbon_cost(brfn, tech, extra)
                 LUC(plot, feed, crop)
                 cultivation_GHG(plot, feed, crop)
                 avg_feed_GHG(plot, feed)
                 site_biomass(brfn)
                 fff(plot, feed);

$onend

*Solve USDA using MIP maximizing profit;

*quant = 0;

*$ontext
Loop run do

*set new demand constraint
fuel_price_base = fprice(run);
fuel_price(term) = fuel_price_base;
*+ sum(state, state_terminals(term, state)*state_price_delta(state));
naptha_price = fuel_price_base*.86;
*pbcost(plot, brfn, feed) = pbcost_base(plot, brfn, feed) + (fuel_price - 3.55*gge_conversion('fahc'))*diesel_use(plot, brfn, feed)/1000;

*xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.75) = 1;

*quant_last_run = quant;

*put update;
*put "Working on " run.tl /
*putclose;

*Resolve for new demand constraint
Solve USDA using MIP maximizing profit;


*$ontext
*Calculate results for export
brfn_production(brfn, tech, extra) = sum(feed, gge_conversion(tech)*conv_eff(feed, tech)*x1.l(feed, brfn, tech, extra))*10;
consumed_biomass(plot, feed, crop) = biomass_consumed.l(feed, plot, crop)*10000;
fff(plot, feed) = sum(crop, consumed_biomass(plot, feed, crop));
site_biomass(brfn) = sum((feed, tech, extra), x1.l(feed, brfn, tech, extra));

loop (plot, feed)$(fff(plot, feed) gt 10) do
         avg_feed_pc(plot, feed) = sum((crop), wta(plot, feed, crop)*biomass_consumed.l(feed, plot, crop) + marg_cost_land(crop)*adopted_acres.l(plot, feed, crop))/sum((crop), biomass_consumed.l(feed, plot, crop));
         LUC(plot, feed, crop) = LUC_CO2(plot, feed, crop)*adopted_acres.l(plot, feed, crop)*10000;
         cultivation_GHG(plot, feed, crop) = cultivation_co2(plot, feed)*consumed_biomass(plot, feed, crop);
         avg_feed_GHG(plot, feed) = (sum(crop, LUC(plot,feed, crop) + cultivation_GHG(plot, feed, crop)))/sum((crop), consumed_biomass(plot, feed, crop));
endloop;

loop (brfn, tech, extra)$(brfn_production(brfn, tech, extra) gt 1) do
         brfn_quant(brfn, tech, extra) = sum(feed, conv_eff(feed, tech)*x1.l(feed, brfn, tech, extra))*10;
         brfn_consumption(brfn, tech, extra, feed) = x1.l(feed, brfn, tech, extra)*10000;
         brfn_electricity(brfn, tech, extra) = sum(feed, electricity(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
         brfn_naptha(brfn, tech, extra) = sum(feed, naptha(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
         feedstock_cap(brfn, tech, extra) = sum(feed, brfn_consumption(brfn, tech, extra, feed))/365/.9;
         capital_cost(brfn, tech, extra) =  (sum(feed, bc(tech)*x1.l(feed, brfn, tech, extra)) + ac(tech)*xi.l(brfn, tech, extra))*10;
         annual_capital(brfn, tech, extra) = CRF*capital_cost(brfn, tech, extra);
         O_M(brfn, tech, extra) = (sum(feed, bo(tech)*x1.l(feed, brfn, tech, extra)) + ao(tech)*xi.l(brfn, tech, extra))*10;

         feedstock_procurement(brfn, tech, extra) = sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot)$(p2b.l(feed, plot, brfn) gt 0.001), avg_feed_pc(plot, feed)*p2b.l(feed, plot, brfn)))
                                                         *10*sum(feed, x1.l(feed, brfn, tech, extra)) / site_biomass(brfn);
         feedstock_transport(brfn, tech, extra) = sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot), pbcost(plot, brfn)*p2b.l(feed, plot, brfn)))
                                                 *10*sum(feed,x1.l(feed, brfn, tech, extra))/site_biomass(brfn);
         fuel_distribution(brfn, tech, extra) = sum((term), T.l(brfn, term, tech)*(btcost(brfn,term)+ terminal_cost(term, tech)))
                                                 *10*sum(feed,x1.l(feed, brfn, tech, extra))/site_biomass(brfn);
         brfn_annual_cost(brfn, tech, extra) = annual_capital(brfn, tech, extra) + O_M(brfn, tech, extra) + feedstock_procurement(brfn, tech, extra) + feedstock_transport(brfn, tech, extra) + fuel_distribution(brfn, tech, extra);
*         marginal_feedstock(brfn, tech, extra) = smax((plot, feed, crop)$(x1.l(feed, brfn, tech,extra) gt 1 and p2b.l(feed, plot, brfn) gt 0.001 and biomass_consumed.l(feed, plot, crop) gt 0.001), (wta(plot, feed, crop) + pbcost(plot, brfn))/conv_eff(feed, tech)/gge_conversion(tech));
*         mc(brfn, tech, extra) = (annual_capital(brfn, tech, extra) + O_M(brfn, tech, extra))/brfn_production(brfn, tech, extra) + marginal_feedstock(brfn, tech, extra) + smax((term)$(T.l(brfn, term, tech) gt 0.1), (btcost(brfn,term)+ terminal_cost(term, tech)))/gge_conversion(tech);
         max_trans_dist(brfn, tech, extra) = 0;
*smax((plot, feed)$(x1.l(feed, brfn, tech,extra) gt 0.1 and p2b.l(feed, plot, brfn) gt 0.1), feedstock_od_table(plot, brfn, 'dist') );

         brfn_GHG(brfn, tech, extra) = (conversion_co2(tech)*sum(feed,x1.l(feed, brfn, tech, extra)) - CO2_elect(brfn, tech)*sum(feed,electricity(tech, feed)*x1.l(feed, brfn, tech, extra)))*10000;
         brfn_CI(brfn, tech, extra) = (brfn_GHG(brfn,tech,extra) + (sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot)$(p2b.l(feed, plot, brfn) gt 0.001), avg_feed_GHG(plot, feed)*p2b.l(feed, plot, brfn)*10000)) +
                                      sum((feed)$(x1.l(feed, brfn, tech, extra) gt 0.1), sum((plot), distance(plot, brfn)*c02_per_ton_mile*p2b.l(feed, plot, brfn)))*10000)*sum(feed,x1.l(feed, brfn, tech, extra))/site_biomass(brfn))
                                         /(brfn_production(brfn, tech, extra)*129.84) ;
         brfn_carbon_cost(brfn, tech, extra) = brfn_production(brfn, tech, extra)*129.84*carbon_price/1000;
         avg_cost(brfn, tech, extra) = (brfn_annual_cost(brfn, tech, extra)- brfn_electricity(brfn, tech, extra)*elec_price + brfn_carbon_cost(brfn, tech, extra))/(brfn_production(brfn, tech, extra) + 0.86*brfn_naptha(brfn, tech, extra));

endloop;

*fuel_production(tech) = sum((feed, brfn, extra)$(x1.l(feed, brfn, tech, extra) gt 0.05), 10*gge_conversion(tech)*conv_eff(feed, tech, extra)*biogenic_energy(feed)*x1.l(feed, brfn, tech, extra));
feedstock_consumption(feed) = sum((brfn, tech, extra), x1.l(feed, brfn, tech, extra))*10000;







*feedstock_consumption(feed) = sum((brfn, tech, extra), x1.l(feed, brfn, tech, extra));
brfn_count(tech) = sum((brfn, extra), xi.l(brfn, tech, extra));
elect_quant = sum((brfn, tech, extra, feed), electricity(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
naptha_quant = sum((brfn, tech, extra, feed), naptha(tech, feed)*x1.l(feed, brfn, tech, extra))*10;
fuel_by_tech(tech) = sum((brfn, extra), brfn_production(brfn, tech, extra) + 0.86*brfn_naptha(brfn, tech, extra));
total_annual_cost = annual_cost.l*10;
brfn_capital_total(tech) =  sum((feed, brfn, extra), bc(tech)*x1.l(feed, brfn, tech, extra)) + sum((brfn, extra), ac(tech)*xi.l(brfn, tech, extra));

*links(plot, brfn, feed) =  p2b.l(feed, plot, brfn)*10000;
*term_links(brfn, term, tech) = T.l(brfn, term, tech)*10;

crops_displaced(crop) = crop_displaced.l(crop)*10000/1000000;
total_land = sum(crop, crops_displaced(crop));
crop_acres_changed(plot, feed, crop) = adopted_acres.l(plot, feed, crop)*10000;
biomass_production(crop) = sum((plot, feed), consumed_biomass(plot, feed, crop)/1000);
total_GHG = total_carbon.l/100;
quant = sum((brfn, tech, extra), brfn_production(brfn, tech, extra) + 0.86*brfn_naptha(brfn, tech, extra));


*Save results
*Model and solver status output file
put status_%scenario%;
put run.tl, "Model status", USDA.modelstat, "Solver status", USDA.solvestat/;

*Summary of national results
put results_%scenario%;
put      "%scenario%", fuel_price_base, quant, total_land, total_GHG,
         loop (crop) do
         put crops_displaced(crop),
         endloop;
         loop (crop) do
         put biomass_production(crop),
         endloop;
         loop (tech) do
         put fuel_by_tech(tech),
         endloop;
put /;

*Site specific supply curve
put results_%scenario%_brfn;
loop     (brfn, tech, extra)$(brfn_production(brfn, tech, extra) gt 1)         do
put "%scenario%", fuel_price_base, brfn.tl, tech.tl, extra.tl, brfn_quant(brfn, tech, extra), brfn_electricity(brfn, tech, extra), brfn_naptha(brfn, tech, extra), feedstock_cap(brfn, tech, extra),
         loop (feed)     do
         put brfn_consumption(brfn, tech, extra, feed),
         endloop
put capital_cost(brfn, tech, extra), brfn_annual_cost(brfn, tech, extra), annual_capital(brfn, tech, extra), O_M(brfn, tech, extra), feedstock_procurement(brfn, tech, extra), feedstock_transport(brfn, tech, extra),
         fuel_distribution(brfn, tech, extra), brfn_carbon_cost(brfn, tech, extra), avg_cost(brfn, tech, extra), max_trans_dist(brfn, tech, extra), brfn_GHG(brfn, tech, extra), brfn_CI(brfn,tech,extra) /
endloop;

$ontext
*Deliveries between plots and biorefinieries
put results_%scenario%_feedstock_links;
loop (plot, brfn, feed)$(links(plot, brfn, feed) gt 100) do
         put "%scenario%", fuel_price_base, plot.tl, brfn.tl, feed.tl, links(plot, brfn, feed)/
         endloop;

put results_%scenario%_fuel_links;
loop (brfn, term, tech)$(term_links(brfn, term, tech) gt 0.1) do
         put "%scenario%", fuel_price_base, brfn.tl, term.tl, tech.tl, term_links(brfn, term, tech) /
         endloop;
$offtext

put results_%scenario%_pc;
loop (plot, feed, crop)$(biomass_consumed.l(feed, plot, crop) gt 0.000100)  do
         put "%scenario%", fuel_price_base, plot.tl, feed.tl, crop.tl, consumed_biomass(plot, feed, crop), crop_acres_changed(plot, feed, crop), LUC(plot, feed, crop), cultivation_GHG(plot, feed,crop)/
         endloop;

if USDA.solvestat = 1 then
         xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.5) = 1;
endif;
*xi.up(brfn, tech, 'extra2')$(xi.l(brfn, tech, 'extra1') gt 0.5) =1;

*xi.lo(brfn, tech, extra)$(xi.l(brfn, tech, extra) gt 0.1) = xi.l(brfn, tech, extra);
*$offtext
endloop;

